Direct observation of hydro dynamic instabilities in driven non-uniform colloidal 

dispersions 
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A Rayleigh- Taylor- like instability of a dense colloidal layer under gravity in a capillary of microflu- 
idic dimensions is considered. We access all relevant lengthscales with particle- level microscopy and 
computer simulations which incorporate long-range hydrodynamic interactions between the parti- 
cles. By tuning the gravitational driving force, we reveal a mechanism whose growth is connected 
to the fluctuations of specific wavelengths, non-linear pattern formation and subsequent diffusion- 
dominated relaxation. Our linear stability theory captures the initial regime and thus predicts 
mixing conditions, with important implications for fields ranging from biology to nanotechnology. 

PACS numbers: 82.70.Dd, 47.57.ef, 61.20.Ja, 05.40.Jc 



Particulate dispersions have long been subjected to ex- 
ternal fields as a means to separate different constituents; 
in particular, sedimentation is important not only for an- 
alytical but also for preparative purposes [1]. For bulk 
systems, successful separation depends crucially upon 
avoiding hydrodynamic instabilities. The development 
of microfluidics [2] has made it possible to exploit the 
suppression of turbulence at small lengthscales in order 
to design novel separation devices [3]; on the other hand, 
this significantly increased stability against mechanical 
perturbations severly limits mixing needed for many 'lab- 
on-a-chip' applications. Often strong external fields [4] 
or complex fabrication [5] are required to produce hydro- 
dynamic instabilities required for efficient mixing. 

Experiments 0] and computer simulations [8] which 
study velocity fluctuations have played a crucial role in 
our understanding of how dispersions respond to exter- 
nal driving fields, in particular to gravity. The motion 
of a solute particle is characterised by a Peclet num- 
ber Pe = td/ts^ which is the ratio between the time 
td it takes a particle to diffuse its own radius and the 
time Ts it takes to sediment the same distance. A Peclet 
number of order unity is the dividing line between col- 
loidal (Pe ^ 1) and granular systems (Pe ^ 1), i.e. Pe 
measures the importance of Brownian motion. All at- 
tempts at a quantitative description of sedimentation to 
date considered a homogenously distributed dispersion 
as the initial state. For preparative purposes, on the 
other hand, starting with a particle-rich layer on top of 
pure solvent is more relevant as it enables the separation 
of particles depending on their sedimentation coefficient. 
However, this configuration is unstable with respect to 



gravity. The particle velocities become correlated, which 
leads to emergent density fluctuations and consequently 
more rapid sedimentation than Stokes' flow alone. It is 
well known that many practical particle concentrations 
develop this Rayleigh- Taylor (RT) like instability. This 
provides an avenue by which the system may be success- 
fully mixed on the one hand, conversely this very mixing, 
leads, chaotically, to a scenario in which separation does 
not occur. For stable separation, it is essential to avoid 
the RT instability. It is possible to use a density gradient 
to counteract the instabilities [ij. 

The 'original' Rayleigh- Taylor instability, which occurs 
if a heavy, immiscible fluid layer is placed on top of a 
lighter one has been intensively studied for the case of 
a simple Newtonian fluid both by theory [9], simulation 
[Tot and experiment, and is observed in granular mat- 
ter [uJ, ll2Lin surface-tension dominated colloid-polymer 
mixtures [l3[ and in a suspension of dielectric particles 
exposed to an ac electric field gradient [14]. 

Here we consider a suspension of colloidal hard spheres 
(without surface tension) of microfluidic dimensions, in 
which we have access to all relevant length scales, from 
the single particle level to the full system. A system- 
atic study of sedimentation in an inhomogenous system 
is presented. We employ three approaches: experiment, 
computer simulation and theory. The experimental real- 
isation is provided by confocal microscopy at the single- 
particle level nisi , while the simulation is a particle-based 
mesoscale technique [Tg*] which captures the direct in- 
teractions between the colloidal particles, and, crucially, 
the solvent which mediates the hydrodynamic interac- 
tions and whose backflow drives the RT instability. Our 
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results at short times are modelled with a linear stability 
analysis 

The RT instability is thought of as a fluctuation in the 
interface between two fluids. Since in a hard-sphere sus- 
pension there is no phase separation, we consider a con- 
tinuous density profile, albeit rapidly varying. To capture 
the lateral fluctuations, we consider the stability of this 
density and associated pressure profile against fluctua- 
tions of wavelength A in a horizontal plane perpendicular 
to gravity. We consider a slit geometry of height L which 
is sketched in Fig. [1^. In the absence of surface tension, 
the fluctuations of all wavelengths are in principle un- 
stable, but short wavelength fluctuations are washed out 
by diffusion of the colloidal particles and so do not grow 
exponentially [itI. fisj. 

Our linear stability analysis, which reveals the stable 
and fast-growing wavelengths of fluctuations, is based 
on a continuum hydrodynamics approach where the col- 
loidal dispersion is considered as an incompressible one- 
component fluid with inhomogeneous mass density p{x) 
and corresponding kinematic viscosity v{x) as obtained 
from the Saito representation pj|. The spatially varying 
density profile is given by p{x) = (j){x)pc + (1 — (j){x))ps^ 
where pc and ps are the mass densities of the colloidal 
particles and the solvent. The colloidal packing frac- 
tion profile (j){x) is an input from an equilibrated sim- 
ulation for inverted gravity. The stability of the initial 
density p{x) and pressure p{x) profiles against pertuba- 
tions dp (X 5p (X ex.-p {i{kyy -\- kzz) -\- n{k)t) with wave 
number k = {ky-\- k'lY^'^ in the yz plane and growth rate 
n is calculated via the linearized Navier-Stokes equations 
resulting in the eigenvalue problem 



n{(p<)' - pPu:,} = {v{u'^' 
-P{f-p'u^ + 2iy'u', 



(1) 



with the spatial derivative . . / = d. . . /dx, the strength 
of the gravitational field g and the fluid velocity field 
in gravity direction Ux{x). For a system confined be- 
tween two rigid walls we impose Ux = along with the 
no-slip boundary conditions dux/dx = at x = 0,L. 
We account for colloid diffusion by the correction term 
n*(/c) = n(k) - Dp [13, [3, with diffusion constant 
D = /cbT/Sttt^sCt (cr colloid diameter) and dynamic sol- 
vent viscosity rjg. 

In our computer simulation, which includes solvent- 
mediated momentum transfer between the colloidal parti- 
cles, we consider a suspension of TV = 15, 048 hard sphere 
particles of mass M immersed in a bath of typically 
Ns = 14, 274, 843 solvent particles of mass m and number 
density = Ns/V. The solvent particles are subjected 



to multi-particle collision dynamics [16|, |20| , which con- 
sists of two steps. In the streaming step, solvent particles 
move ballistically for time St. In the collision step, parti- 
cles are sorted in cubic cells of size a, and their velocities 
relative to the center-of-mass velocity of each cell are ro- 
tated by an angle a around a random axis. We employed 
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FIG. 1: a, A schematic illustrating the spatial parameters a, A 
and L. b-e, Simulation snapshots of a system which contains 
N = 33, 858 colloidal particles and Ns = 32, 118, 397 solvent 
particles (not displayed) in a simulation box with dimensions 
L/a = 18 and Ly/a = Lz/a = 81. The value of the Peclet 
number is Pe = 1.6. b-d, Time series of the system at time 
t/rs = 3.2 (b), 6.4 (c), 9.6 (d). The snapshots are slices 
of thickness 2a done in the xy plane, e, Slice of thickness 
2a in the yz plane at time t/rs = 9.6. The height of the 
yz plane is x/L = 2/3, as indicated by the dashed line in 
(d). f-i, Experimental realisation of the Rayleigh- Taylor- like 
instability, f-h, Time series of images taken with a confocal 
microscope in the xy plane for the parameters = 0.15, Pe = 
1.1 and L^/cr = 18 at times t/rs = 1.43 (f), 5.5 (g), 11.22 
(h). i, Slice in the yz plane at a height x/L = 2/3 (indicated 
by the dashed line in (h)) at time t/rs — 11.22. In (f-i) the 
scale bar denote 40 /xm. 



the parameters 6t = ^.2^ma^ /ksT ^ a = 37r/4, n^a^ = 5 
and M = 167m in order to achieve the hierarchy of time 
scales and the same hydrodynamic numbers as in the 
experiment, see Ref. [l, [2l| for details. To enforce no- 
slip boundary condition on the colloid surface and the 
confining walls a stochastic-reflection method [22[ is ap- 
plied. Statistical averages for time-dependent quantities 
are performed over 200 independent configurations. 

In our single-particle level confocal microscopy experi- 
ments we used polymethylmethacrylate colloids sterically 
stabilised with polyhydroxy-stearic acid. The colloids 
were labeled with the fluorescent dye coumarine and had 
a diameter a = 2.8 pm with around 4% polydispersity as 
determined by static light scattering. To almost match 
the colloid refractive index we used a solvent mixture 
of cis-decalin and cyclohexyl bromide (CHB), which we 
tuned to yield different Peclet numbers, owing to changes 
in the degree of density mismatch between colloids and 
solvent. Specifically, Pe = 1.1 and Pe = 2.4 correspond 
to 80% and 87.5% CHB by weight respectively. The char- 
acteristic time to diffuse a radius td ^ 29 s. The data 
were collected on a Leica SP5 confocal microscope, fitted 
with a resonant scanner, at a typical scan-rate of around 
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10 s per 3D data set. Prior equilibration was achieved by 
placing the suspension overnight such that it sedimented 
across a thin (typically 50 jam) capillary. The capillary 
was then inverted, and the evolution under sedimentation 
was followed. 

We begin our discussion by presenting snapshots of 
the system, in Fig. [T}3-e from computer simulation, and 
in Fig. [Tf-i from confocal-microscopy. The similarity is 
remarkable, and we note that, at the very least, our sim- 
ulation qualitatively reproduces the experiment. For a 
quantitative comparison, we consider the dispersion rela- 
tion of wavenumber against growth rate in Fig. [2)3. The 
time evolution in the development of the RT instability 
with a characteristic wavelength is clear. While snap- 
shots in the gravity plane (Fig. [TJd, c, d, f, g, and h) 
illustrate the overall process of sedimentation, snapshots 
in the horizontal yz plane show the transient pattern or 
network-like structure that results from the RT instabil- 
ity (Fig. [1^ and i). At later times, the network structure 
decays and a laterally homogenous density profile devel- 
ops where the colloids start to form a layer at the bottom 
of the cell which becomes more compact with time. The 
time evolution is shown in detail in the Movies 1-4, see 
EPAPS Document No. []. 

The linear stability analysis predicts the existence of 
the initially fastest growing wavelengths in the RT insta- 
bility. We plot the results of the linear stability analysis 
for a range of slit widths L keeping Pe fixed, and for a 
variety of Peclet numbers keeping L fixed in Fig. [2^1, b 
and c respectively. The dimensionless growth rates, utd^ 
are plotted as a function of wave number kcr = 27r(j/A, 
where A is the wavelength of the fluctuations as indicated 
in Fig. [1^. Without diffusion, fluctuations at all wave 
numbers are unstable, as shown by the solid lines in Fig. 
[2^, b and c. Due to diffusion, we find that growth rates 
at higher wavevectors are suppressed as expected, i.e., 
diffusion destroys the Ray leigh- Taylor instability at suf- 
ficiently small wavelengths. We find excellent agreement 
between the theory with diffusion and both simulation 
and experimental data, up to /ccr ~ 1, which is surpris- 
ing for a coarse-grained continuum description. With 
decreasing wall separation L, the growth rate Umax de- 
creases and kmax increases, see Fig. [2^. Since the fluid 
velocity in the gravity direction decreases as e~^^, where 
X is the distance from the interface, only lon^ wavelength 
undulations feel the presence of the walls [9, 23]. Figure 
[2)3 and c show that for fixed L, driving the sedimentation 
more strongly by increasing the Peclet number leads to 
an increase in the wave number of the fastest-growing un- 
dulation kmax^ and the corresponding growth rate rimax- 

So far we have considered only the linear regime of the 
instability, which is valid at small times, when the ampli- 
tude of the fluctuations is smaller than the wavelength. 
Our experiments and simulations permit detailed access 
to all relevant time- and length-scales in the non-linear 
regime, where the colloids form foam-like structures in 
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FIG. 2: a-c Growth rate utd versus wave number ka = 
27va/X. a, Simulation results of n{k) for different wall sep- 
aration distances L/a — 18, 12, 9 and fixed Pe — 1.6. b, Sim- 
ulation and experimental results of n{k) for different Peclet 
numbers Pe = 4.8,1.6,1.1,0.8,0.4 and fixed L/a = 18. c, 
Experimental results of n(k) for different Peclet numbers 
Pe — 2.42, 1.21 and fixed L/a — 36. The open symbols are 
the results obtained from simulation, filled symbols are exper- 
imental results, solid lines represent the solutions from the in- 
stabilty analysis and the dashed lines are the same numerical 
solutions plus the diffusion correction. First moment of the 
colloid density (x) / L (d), (f), (h) and second moment of the 
colloid density ax/L (e), (g), (i) versus time t/rs. Solid lines 
indicate simulation data whereas the dashed lines indicate 
experimental data. d,e, L/a = 18,12,9 and Pe = 1.6. f,g, 
Pe = 4.8, 1.6, 1.1, 0.8, 0.4 and L/a = 18. h,i, Pe = 2.42, 1.21 
and L/a — 36. 



the (confined) xz plane (Fig. [Tt,g) and a network-like 
structure in the yz plane (Fig. [l^,i) appears. Appar- 
ently, both continue to exhibit the characteristic length 
scale \max = '^Tr/kmax of the fastest growing wavelength 
in the linear regime. 

In order to quantify the different regimes of the in- 
stability we use the first moment of the density (x), i.e. 
the centre of mass of the colloid coordinates and the sec- 
ond moment of the density = (x^) — {x^ . Here, (x) 
is a measure of the degree of sedimentation, while ax 
quantifies the extent to which the instability spreads out 
the colloids in the gravity direction. Three regimes are 
clearly visible in Fig. [2]d, f, h, e, g and i. Firstly we 
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FIG. 3: Time series of images taken with a confocal micro- 
scope in the xy plane for the parameters = 0.15, Pe = 2.42 
and L^/cr = 36 at times t/rs = 3.1 (a), 8.1 (c), 24.6 (d) 
and 56 (e). The crystahine layers are clearly visible, (b) is 
a slice in the yz plane approximately in the middle of the 
colloidal crystal in (a). The secondary instability occurs at 
t/rs ~ 20 — 30, see (d). In (a,c-e) the scale bar denotes 40 
/im and 20 /xm in (b). 

find the linear regime in which the flat interface devel- 
ops undulations and hence {x) slowly decreases and cjx 
slowly increases, secondly the non-linear regime where 
'droplets' of colloid-rich material fall to the bottom and 
therefore {x) sharply decreases and (Jx sharply increases, 
and thirdly the regime in which the colloids start to form 
a layer at the bottom of the cell which becomes more 
compact with time under settling as can be seen from 
the slow decrease of both {x) and (Jx- Clear agreement 
between simulation and experimental data can be seen 
from Fig. [2] f and g. 

In the case of a rather large slit width L/a = 36, there 
is a sufficiently high sediment for a region of colloidal 
crystal to form, see Fig [3] and EPAPS Document No. [] 
for Movie 5. Since the crystal has a finite (albeit small) 
yield stress, the only flow we observe initially occurs in a 
thin fluid layer between the crystal and the lower solvent 
region via narrow vertical tubes, in marked contrast to 
Fig. [T]b-d. The crystal melts layer by layer until flnally 
it becomes sufficiently thin that it peels off the wall in a 
second instability, which leads to a change of slope for the 
Pe = 2.42 line in Fig. Eh and i at t/rs ~ 20-30, see Fig. 
[3] d, until most of the particles have sedimented down 
(Fig. [3]e). This observation of driven surface melting 
at the single particle level has the potential to provide 
new insight into this poorly understood phase transition 
under non-equilibrium conditions. 

Using state-of-the-art simulation and experimental 
techniques, we have presented a quantitative analysis 
of a hydrodynamic instability in a colloidal system at 



a microfluidic lengthscale. Our results show excellent 
agreement between experiment and simulation, showing 
that the latter accurately describes the fundamentally 
and practically important phenomena caused by hydro- 
dynamic instabilities. Furthermore, by employing a sim- 
ple theoretical treatment to the initial linear behaviour, 
we flnd considerable predictive power. The theory can 
flexibly be used to predict conditions for separation and 
mixing. We also note that the theory reveals even the 
length scale of the network structure that results from 
the instability. We flnally emphasise that the full access 
and accuracy to all relevant length scales in this prob- 
lem allowed for the observation of novel phenomena, not 
yet explored further, such as the inverse gravity induced 
crystal melting. 
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